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Abstract 

Chronic diseases impose a tremendous global health problem of the 21st century. Epidemiological and public health 
models help to gain insight into the distribution and burden of chronic diseases. IVloreover, the models may help to plan 
appropriate interventions against risk factors. To provide accurate results, models often need to take into account three 
different time-scales: calendar time, age, and duration since the onset of the disease. Incidence and mortality often change 
with age and calendar time. In many diseases such as, for example, diabetes and dementia, the mortality of the diseased 
persons additionally depends on the duration of the disease. The aim of this work is to describe an algorithm and a flexible 
software framework for the simulation of populations moving in an illness-death model that describes the epidemiology of 
a chronic disease in the face of the different times-scales. We set up a discrete event simulation in continuous time involving 
competing risks using the freely available statistical software R. Relevant events are birth, the onset (or diagnosis) of the 
disease and death with or without the disease. The Lexis diagram keeps track of the different time-scales. Input data are 
birth rates, incidence and mortality rates, which can be given as numerical values on a grid. The algorithm manages the 
complex interplay between the rates and the different time-scales. As a result, for each subject in the simulated population, 
the algorithm provides the calendar time of birth, the age of onset of the disease (if the subject contracts the disease) and 
the age at death. By this means, the impact of interventions may be estimated and compared. 
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Introduction 

Chronic diseases impose a tremendous global health problem of 
the 21st century. The World Health Organization estimates that 
63% of all deaths in 2008 were caused by chronic diseases [1]. 
Besides taking measures in politics and society, research efforts are 
needed to oppose this threat. In studying the characteristics of 
chronic diseases from a public health perspective, it is often 
important to consider different time-scales [2]. The age of the 
subjects in a population is a risk factor for many diseases. Changes 
in life-style and medical care influence the risk of contracting and 
dying from the disease on a secular scale. Thus, the incidence rate 
of the chronic disease as well as the mortality rates (with and 
without the disease) depend on the calendar time. Moreover, the 
mortahty of people with the disease often depends on the duration 
since its onset. Examples are diabetes [3,4], dementia [5], 
depression [6], and systemic lupus erythematosus [7]. For 
decision-makers all time-scales may be important. 

As a hypothetical example, consider the question which health 
programme to choose from two possibilities A and B if the 
outcome of interest is the gain of life-years. Possibility A is known 
to decrease the incidence of the disease by 15%, and possibility B 
lowers the mortality of those having the disease for more than ten 
years by 50%. The decision depends on several factors. If, on the 



one hand, the incidence rate is non-zero for children only and the 
birth rate in the population is low, possibility A may have little 
effect with respect to the gain of life-years. If, on the other hand, 
the chronic disease has very few people reaching ten-year survival 
after onset, programme B can be nearly useless. In problems 
similar to the example, the decision-maker may face a complex 
interplay of epidemiological and demographical considerations. 

The aim of this work is to describe an algorithm and a flexible 
software framework for simulation of populations moving in a 
multi-state model (illness-death model) that describes a chronic 
disease. The simulation takes into account the different time- 
scales: calendar time, age, and duration of the disease. Although 
simulations using multi-state models are subject to recent 
textbooks [8], to our knowledge no algorithm has been described 
that incorporates the effects of all the different time-scales. 

Methods 

A popular framework for studying irreversible diseases is the 
iUness-death model (IDM) consisting of the three states Normal, 
Disease and Death, [9-11]. The associated transition rates, 
synonymously densities (in units "per person-time", not to be 
confused with risks or probabilities [12]), are the incidence /, and 
the mortality rates OTq and m\ (Figure 1). In general, these rates 
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Figure 1. Three states model of normal (healthy), diseased and 
dead subjects. The transition rates may depend on calender time t, 
age a, and in case of mi also on the duration d of the disease. 
doi:1 0.1 371/journal.pone.01 06043.g001 

depend on calendar time (t), age (a) and in case of in[ on tlie 
duration of the disease (d). 

This article presents a method for simulating populations 
moving in the IDM. The motivation for the algorithms comes 
from analytical epidemiology where relations between common 
epidemiological measures are studied. Examples for those 
measures are the prevalence, the duration of a disease, the age 
of onset (or diagnosis), and lost life-years (due to the disease). A 
typical question may be: what is the mean age of diagnosis of 
subjects born in a certain time period? What is their mean age at 
death? Another interesting aim is the estimation of the incidence 
rate / from cross-sectional information. At a specific point in time 
t', each of the subjects j=l,...,n, has a unique "status". 
Neglecting those who are unborn or dead at t', the status in the 
IDM is either normal (non-diseased) or diseased. Thus, the status 
can be seen as a binary random variable, and data of this kind are 
typically called current status data [13]. The current status is 
closely linked with the incidence / and the mortalities ma and m\ 
before t' . Estimating the incidence from current status data, for 
example, has been a topic in research for decades [14]. The 
framework presented here may be useful in this field. 

Overview of the simulation algorithm 

The simulation is a microsimulation, i.e., it treats each person in 
the population as an autonomous unit. For each person, indexed 
7,7=1,... ,M, the relevant events diagnosis and death are 
simulated. This is accomplished in two steps: 

1 . Contracting the disease or dying without the disease is 
modelled as competing risk [11]. Given the time t^ of birth 
of person 7, the cumulative distribution function Fj^ ' of the /mi 
failure time T^^ is 



i^W(r) = l-exp 



-T,T Idr 



The Xtrm first failure time T'f'* refers to the time of diagnosis or 
death without disease and is measured in time units after birth 
of person 7. Thus, T^^ is the age at which the first transition 
from the state Normal occurs. Given that a transition occurs at 
for person 7, then the odds of moving into state Disease 
versus moving into state Death is 



1 



2. If the event at T^~^ is the death (without the disease), the 
simulation for person 7 is finished. If, however, the event is the 
diagnosis of the disease, the "second failure time" 7"^' to death 



The next section describes in detail how the integrals in 
Equations (1) and (2) may be calculated in the simulation. After 
calculating the integrals, the question arises how the times T\ and 
Ti can be obtained from F\ and Fj. This is done by the inverse 
transform sampling method: Let F be a cumulative distribution 
function and Me(0,l). For F-\u) : =inf {x\F{x)>u} it holds: If 
(7 is a uniform random variable on (0,1), then F^^(U) follows the 
distribution F. Thus, the simulation of Ti and Tj is easy, if a 
random number generator for U such as runif in R is available. 

For each of the n persons in the population we store four pieces 
of data: 

1 . a unique identifier 7, 

2. the date fp ' of birth (dob) of person 7, 

3. the age at diagnosis (adi) of person 7, and 

4. the age at death (ade) of person 7. 

If the person 7 does not contract the disease, the age at diagnosis 
adi is set to NA (missing). In summary, we get the Algorithm 1. 
Algorithm 1 Simulation of populations moving in the IDM 
for 7=1 to n do 



2 
3 
4 

5: 
6: 
7: 
8: 
9: 
10: 
11: 
12: 
13: 



dob 



calculate event time T^^ according to Equation (1) 

simulate type of event that has happened at 7""' by 
Equation (1) 

if event is diagnosis then 



r(/) 



adi 



r(j) 



calculate time 7"^ of death using Equation (2) 



ade <- 
else 
adi <- 
ade <- 
end if 
write j , 
14: end for 



r(J}_ 



NA 



</") 



dob, adi, ade to file 



Calculating line integrals 

In this section the calculation of the integrals in Equations (1) 
and (2) is described. The situation in which analytical expressions 
for these integrals exist, is straightforward. The first simulation in 
the next section is an example. However, in real world 
applications, analytical expressions for the integrals are rarely 
given. For convenience, mathematical functions (e.g., splines) may 
be fitted to the data and integration is accomphshed with the fitted 
functions. Since the aim of this work is a flexible way of treating 
the incidence and mortality rates, we assume that the rates are 
given as numerical values on a regular grid. Here we focus on the 
most general case, which is characterized by: 

1 . None of the time-scales t,a, and d is negligible, and 

2. the values of i,mQ,m\ are given as data points only. 

For many chronic diseases, we think this is the most relevant 
case: The mortality rates mo and m\ depend on a and t. Since age 
is a risk factor for many diseases, the dependency on age is 
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obvious. Healthier life-style and medical progress in many 
countries lead to secular trends in mo and ini . In addition, disease 
duration d is likely to have an impact on mi in many chronic 
diseases. Thus, none of the time-scales is negligible, which is 
covered in the second and third example. 

In the most general case, the integrands i,mo, and mi are given 
by data points only. We assume that the numerical values are 
located on a regular grid. The grid is two-dimensional in case of i 
and mo, which depend on two time-scales t and a; and the grid is 
three-dimensional in case of mi which depends on t,a, and d. 

In event history analysis [2], a useful concept is the Lexis 
diagram, which is a co-ordinate system with axes calendar time t 
(abscissa) and age a (ordinate). The i-dimension sometimes is 
referred to as period. Each subject is represented by a line segment 
from time and age at entry to time and age at exit. Entry and exit 
may be birth and death, respectively, or entry and exit in a 
epidemiological study or clinical trial. There are excellent and 
extensive introductions about the theory of Lexis diagrams (see for 
example [9,15,16] and references therein), which allows us to be 
short here. In irreversible diseases, the common two-dimensional 
Lexis diagram with axes in t- and a-direction may be generalized 
to a three-dimensional co-ordinate system with disease duration d 
represented by the applicate (i.e., the z-axis). If a subject does not 
contract the disease during lifetime, the life line remains in the t-a- 
plane parallel to the line bisecting abscissa and ordinate. In other 
words, the life line for the time without disease is parallel to 
ei : =(1,1,0) (where the triple (t,a,d) denotes the co-ordinates in 
time, age and duration direction, respectively). However if at a 
certain point in time E the disease is diagnosed, the life line 
changes its direction, henceforth runs parallel to 62 : =(1,1,1). 

The situation is illustrated in Figure 2. The life lines of two 
subjects are shown in the three-dimensional Lexis space. At birth 
(denoted Br,v= 1,2) both subjects are disease-free; both life lines 
are parallel to 6]. The first subject contracts the disease at E. 
Henceforth, the life line is parallel to £2 until death at Di. The 
second subject remains disease-free until death at D2. 

Having the concept of the Lexis diagram at hand, we observe 
that El and Ej in Equations (1) and (2) are line integrals in the 



Lexis space. We start with calculating the first failure times Ti . For 
subject j the associated life line starts at (t,d) = (?o'>0)- We chose an 
age CO > 0, when it is certain that a transition to one of the states 
Disease or Death has occurred, say (u= 150 (years). For calculating 



E'f\ we trace the hypothetical life line from Bj 



--{t^\0) to 



Dj : =(?Q +cu,co). Thus, the hypothetical life line has a repre- 
sentation 

Cj : Bj + aiiDj-Bj),oie[Q,l]. 



Following the life line is related to the method of ray-tracing in 
the field of computer graphics, where efficient algorithms for this 
purpose exist. In Siddon's algorithm [17], the key idea is to follow 
Cj by calculating intersections with volume elements (voxels), 
which form a regular partition of the Lexis space. Let 

^* = {o:("(p)|p=l,...,pW} 

with 0 = a'''(l)< ... <0£^'(P^')=1 be a parametrization of the 
points where £j intersects the voxel faces plus the start and end 
points Bj and Dj. Details for the calculation of are described in 
the supporting information to this article. The parametrization A* 
is ideally suited for approximating the integral in Equation (1) by 
the trapezoidal rule [18]. The reason lies in the fact that in 
calculating Ff '(co) the values i^j'' (^?{(' -|-aW(p)co) ,p= 1, . . . ,pW, 

are a by-product. Algorithm 2 shows the necessary steps. 
Algorithm 2 Calculating Ei 
for 7 =1 to n do 



2 
3 
4 

5 
6: 



calculate ^; = {o(('>(p)|P= 1, • • ■ ,^^'} 
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Figure 2. Three-dimensional Lexis diagram with two life lines. Abscissa, ordinate and applicate (z-axis) represent calendar time t, age a and 
duration d, respectively. The life lines start at birth 5,. and end at death D,.,y = 1,2. The first subject (blue line segments) contracts the disease at E. 
Then, the life line changes its direction. The second subject (red line segment) does not contract the disease, the life line remains in the f-a-plane. 
doi:1 0.1 371/journal.pone.01 06043.g002 
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7: for;7 = 2 to do 

8: Zp^oi^Hpyoj 

9: fp^i{t'i' + 

1 



10: ep^ep-i+-iTj, 
'(J) fJJ) _ 



^p-&(fp+fp-l) 

12: end for 
13: end for 

Since the values of / and mo are given on the voxel grid only, the 
calculation oifp,p=l,... ,P^\ needs bilinear interpolation of the 
values of the adjacent voxels [19]. 

After preparing F^^ j= \ , . . . ,n, the values of the times T^^ can 
be calculated by the inverse transform sampling method. Since we 

have fJ'' calculated at points C/j : = ?o' + ■^/"/'= 1' ■•• the 
inverse transform sampling would yield only those C^. A better 
accuracy can be obtained by (linear) inverse interpolation [18]. For 

Cp-Cp-i 



te(C^_i,y,p = 2,...,pW, let^ 



Then, it holds 



(i-a-^-["(c,-,)+c-^-r(g 



Thus, for Me(0,l) drawn from a uniform distribution, we find 
the unique pe{2, . . . ,P^^, such that Fl'\Cp_i)<u<F\'\Cp), and 
set 



as the inverse {^i ^^ (")• 

For those subjects / who contract the disease, the associated 
F^\'\T^ ') can be derived in a similar way as in Algorithm 2. The 
associated line segment starts at (/,a,^/) = (7q' + rj' *,?'}'', 0). 
Again, a hypothetical maximal disease duration oj' is assumed, 
say to' = 80 (years), such that the line segment ends at 
(?,a,rf) = (?Q ' + T'l' ' + oj',7'j' ' +ft)',co'). Thus, the line segment is 
parallel to e2 = (l5l,l)- The Siddon algorithm computes the 
corresponding set of intersections with the voxel grid accordingly. 



The ages ' at death with disease are obtained from Algorithm 2 
mutatis mutandis. The interpolation of Wi needs to be trUinear. 

Examples 

This section shows the results of different simulation settings. 
The associated R [20] code is provided with this article. The first 
simulation is about a hypothetical chronic disease with rates 
i,mQ,m\ only depending on a. The corresponding age-specific 
prevalence can be calculated analytically, which allows cross- 
checking the results of the simulation. In the second example, 
another hypothetical disease is treated with mortality rates 
depending on t,a, and d. Here we use the ray-tracing approach 
in the Lexis diagram. Again, the outcomes of the simulation are 
compared with analytical results. The third simulation is about a 
relatively rare rheumatic disease. A hypothetical birth cohort of 
1 00000 women is followed from birth to death and examined if the 
disease is diagnosed. Finally, the last simulation demonstrates 
applicability in the context of medical decision-making. Two 
interventions are compared with respect to the outcome life-years 



Simulation 1: Analytical example 

In the first simulation, only one time-scale is involved. Assuming 
»io(a) = »?i(fl) = exp(— 10 + O.lfl), and 



i{a)-- 



0 

0.005 
1-0.005 (fl-20) 



for a < 20 
for a > 20, 



one can show that the age-specific prevalence p(a) is given by 
/>(a) = 0.005(a — 20)^ [21]. The notation x+ means the positive 
part of X : x+ =max(0,x). The integrals of /(fl),mo(a), and m\{d) 
can easily be expressed analytically. Figure 3 shows the age- 
specific prevalence in a simulated cohort with « = 200000 subjects. 
For comparison, the true (analytical) prevalence is depicted as a 
solid line. The result of the simulation agrees very well with the 
analytical age profile of the prevalence, which indicates the 
correctness of the implemented R code. 

Simulation 2: Mortality depends on duration 

The second simulation is about a hypothetical disease with all 
time-scales t,a and d playing a role. We aim at calculating 
epidemiological measures that describe the population of the 
diseased. For example, age of onset, mean duration of the disease 
and age at death may be important to plan resource allocation 
(e.g., inpatient facilities). 

In each of sixty consecutive years 7 = 0, . . . ,59, 2000 persons are 
born and followed from birth to death. The incidence of a 

(fl-30) + 



hypothetical chronic disease is assumed to be i(t,a)- 



3000 



the age-specific mortality rate of the non-diseased is chosen to be 
Wo{7,a) = exp(- 10.7 + 0.1fl + rln(0. 998)) and the mortality of the 
diseased is m\{t,a,d) = mQ{t,a)-{QM(d -if + \). In total, 42299 
of the 120000 simulated persons contract the disease. The 
simulated data allows the derivation of important epidemiological 
measures. For example, the histograms of the age of onset and age 
at death are shown in Figure 4. 
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Figure 3. Theoretical and simulated age-specific prevalence. 

Simulation 1 comprises k = 200000 persons. The resulting age-specific 
prevalence (black crosses) is compared to the analytically calculated 
prevalence (blue solid line). The example shows the very good 
agreement between the simulation and the theoretical results. 
doi:1 0.1 371/journal.pone.01 06043.g003 
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Figure 4. Histograms of the age of onset and age at death in a hypothetical chronic disease. In Simulation 2 the empirical distributions of 

the age of onset (left) and the age at death of the diseased persons (right) are estimated. 

doi:10.1371/journal.pone.0106043.g004 



The median age at death of those who contracted the disease is 
77.20 (years) whereas the median age at death of those without the 
disease is 79.67 (years). The median duration of the disease in the 
42299 ill subjects is 12.69, the interquartile range is 7.46-17.60 
(years). 

Finally, we can cross-check the results of the simulation by 
comparison with an analytical calculation. In year / = 1 00, exacdy 
76548 persons are alive, 8802 of those having the hypothetical 
disease. Figure 5 shows the age-specific prevalence at ?= 100. The 
black lines indicate the prevalence of several age groups together 
with 95% confidence bounds as given by the simulation. The blue 
line represents the prevalence calculated analytically by the exact 
formula in ([9], Section 7.2). The results agree quite well within 
the confidence bounds. 

Simulation 3: Systemic lupus erythematosus in the UK 

The third simulation is about a hypothetical cohort of 100000 
women born in the United Kingdom in 1945. The disease under 
consideration is systemic lupus erythematosus (SLE), which is an 
autoimmune disorder with significant morbidity and increased 
mortality. Women are more often affected than men; in terms of 
age-standardized incidence, the ratio is 5:1 [22]. The peak of the 
age-specific incidence in British women is located at the age of 
about 50. Since Somers et al. could not fmd a significant secular 
trend in the age-standardized incidence in 1990-1999, we assume 
that the incidence does not depend on /, [22]. The aim of this 
example is to study the interplay between age of onset and the 
mortahty of the diseased women in a realistic setting. 
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Figure 5. Calculated and simulated prevalence in Simulation 2. 

If we mimic a cross-sectional study at year t= 100, we obtain the age- 
specific prevalence as indicated by the black bars (with 95% confidence 
bounds). For comparison the analytically calculated age-specific 
prevalence is shown as blue line. 
doi:1 0.1 371/journal.pone.01 06043.g005 
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Mortality mo has been taken from the official life tables of the 
British Office for National Statistics [23]. The calendar time trend 
in »?o has been simplified by assuming that the yearly decrease in 
mortality is 0.36% for all age groups. This is the geometric mean 
of the decrease in female mortality taken over all reported age 
groups in the past 60 years. Mortality nii of the women with SLE 
has been modelled according to [7], which takes into account 
several covariables; sex, age, duration of SLE, calendar time of 
diagnosis and geographical region. Unfortunately, an interaction 
analysis of these factors has not been reported, which forces us to 
make assumptions. We assume a multiplicative model of the 
impact of sex, age, region and SLE duration: 

mi(t,a,d) = c-mo(t,ayHiige{a)'U^m(d). 

The constant c reflects the impact of sex and the region, 77age 
and i/dur are hazard ratios. The hazard ratio i/dur of SLE 
duration has been interpolated affine-linearly on a logarithmic 
scale from the published values and backtransformed. The age 
dependency of //age(a) has been treated similarly. Due to the weak 
(and possibly insignificant) effect, the impact of the calendar time 
of diagnosis has been neglected. 

From the 100000 women in the simulation, 513 contract SLE. 
This corresponds to a lifetime risk of about 0.5%. For comparison, 
a recent publication about women in the US estimated 0.9% [24]. 
The study [24] included women with African ancestors, who are 
known to have a higher risk. The median age of onset in the 
simulation is 46.07 years, which is in nearly perfect agreement 
with the result of 46. 1 years based on another simulation using the 
same data but a different method [25]. Table 1 shows the 
interplay between age of diagnosis and duration of SLE in the 
simulation. 

From Table 1 it is apparent that ten-year survival is not easy to 
reach. Substantial loss of life time is also indicated by the age at 
death: Median age at death for the diseased and the non-diseased 
women is 61.4 and 77.1 (years), respectively. This indicates a 
considerable loss of lifetime in the population of the diseased 
compared to the non-diseased. Fortunately, this situation has 
changed in the last years with better medical care for SLE patients. 
Especially the introduction of optimized treatment regimes in the 
last decade lead to an enormous reduction of mortality [26]. These 
effects have not been included in the simulation. 

Simulation 4: Effectiveness of two interventions 

We take Simulation 2 as the basecase and compare the 
effectiveness of two hypothetical intervention programmes A and 
B with respect to the total gain of life-years. Assume that the 
primary prevention programme A gradually lowers the incidence 
rate for all a>30 in the calendar years ?e[0,10] by 15% and 
remains at the 85% level for /> 10. Programme B is assumed to 
reduce the mortality mi of the diseased people gradually by 50% 
starting after six years with the disease. B could be achieved, for 
example, by listing a new drug on the formulary. After running the 
simulation, we find that in total the 120000 simulated persons in 
the basecase, in intervention A, and intervention B, have 9275008, 
9348941, and 9351338 life-years, respectively. Hence, intervention 
A and B yield about 73900 and 76300 life-years more than the 
basecase. With respect to the chosen outcome parameter, gain of 
life-years, programme B turns out to be superior to programme A. 
Thus, primary prevention by programme A, in this example, is less 
effective than lowering the mortality of the diseased. This result is 
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Figure 6. Esimtated densities of tKie gain in life-years in two 
intervention programmes. Simulation 4 compares two intervention 
programmes in a population of 120000 persons. The estimated 
densities of the total gain in life-years (compared to a basecase) are 
shown as red and blue curve for programmes A and B, respectively. 
doi:1 0.1 371 /journal.pone.01 06043.g006 

hardly predictable without a simulation and demonstrates the 
usefulness of our algorithm in decision making. 

To obtain the empirical distributions of the total gain of life- 
years in both interventions, 5000 random samples of size 120000 
have been drawn (with replacement from the simulated 120000 
persons) [27]. Figure 6 shows the resulting densities and the 
superiority of programme B. 

Conclusions 

This article is about simulating populations in an illness-death 
model consisting of the three states Normal, Disease, and Death. 
The relevant events, diagnosis and death, in the most general case 
depend on three time-scales: calendar time, age and disease 
duration. 

After birth of a subject in the population, two cases may occur: 

1 . the subject dies without the disease, or 

2. the subject contracts the disease and dies with the disease. 

Both situations can be represented in the Lexis space, a 
common tool in event history analysis. In the first case, the life line 
is solely located in the time-age-plane. In the second case, the life 
line changes its direction after onset of the disease, which allows to 
model the duration of the disease. In many diseases, the duration 
plays an important role for the mortality. Beside systemic lupus 
erythematosus as treated above, other diseases such as diabetes 
[3,4], depression [6] and dementia [5] may serve as examples. 

In the most general case, the simulation requires to numerically 
solve line integrals. Therefore, a synthesis between a ray-tracing 
technique and numerical integration is exploited. The method 
provides a fast way to follow the individual life lines in the Lexis 
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